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Abstract. We consider regularization methods of Kaczmarz type in connection 
with the expectation-maximization (EM) algorithm for solving ill-posed equations. For 
noisy data, our methods are stabilized extensions of the well established ordered-subsets 
expectation-maximization iteration (OS-EM). We show monotonicity properties of the 
methods and present a numerical experiment which indicates that the extended OS-EM 
methods we propose are much faster than the standard EM algorithm. 
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1 Introduction 

The expectation-maximization (EM) algorithm provides approximations for maximum likeli- 
hood estimators of problems with incomplete or noisy data, which is the usual framework 
when dealing with inverse or ill-posed problems. In particular, the EM algorithm for Pois- 
son models is well known for its applications to astronomical imaging and to PET (positron 
emission tomography) - see, e.g. [19], [21]. 

In this work we address inverse problems modeled by operator equations which admit non- 
negative solutions, with the aim of approaching them by combined EM-Kaczmarz strategies. 

We begin our study by considering the operator equation 

Ax = y, (1.1) 
where A : L^{il) L^{^) is a Predholm integral operator of the first kind 

{Ax){s) = / a{s,t)x{t)dt, sGS. (1.2) 
Jn 
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Nonnegative solutions of (jl.ip can be determined by finding minimizers of the functional 
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Formally, the first order necessary condition for such a minimizer reads as 
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If the assumption (A* l)(t) = 1 is satisfied, a solution of (|1.3p can be obtained by solving the 
corresponding multiplicative fixed-point equation 
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The fixed-point equation (jl.4|) motivates the definition of the EM algorithm, see [191 Ell El El 
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i.e. an explicit iterative method for solving (jl.4p . 

The OS-EM (ordered subsets - expectation maximization) iteration was introduced in [8] 
as a computationally more efficient alternative to the original EM iteration for the discrete 
case. The main idea is as follows. The data y are grouped into an ordered sequence of subsets 
(or blocks) yj. An iteration of OS-EM consists of a single cycle through all the subsets, in 
each subset updating the current estimate by an application of the EM algorithm in that 
data subset. This strategy can be connected to the Kaczmarz type iterative methods recently 
investigated in [H [H [H [H] for approaching systems of integral equations. 

In order to extend the OS-EM method to infinite dimensional settings, we first group 
the data y into N blocks yj := y\-^^, where Sj C S are not necessarily disjoint and satisfy 
S = So U • • • U Sjv_i. Then equation (jl.ip is decomposed into a system of integral equations 
of the first kind 

A,x = yj, j = 0,...,N-l, (1.6) 

where the Fredholm integral operators Aj : L^{Q) — > L^{T,j) correspond to blocks of A and 
are defined by 



iAjx){s) 



aj{s, t) x{t) dt , 



(1.7) 



with aj := a|nxEj- Notice that x is a solution of (jl.6p if and only if x solves (II. Ij) . 

In order to simplify notation, we drop the indices of the domains and simply write 
Aj : L^{Q) L^i'S) and yj G L^(S). Thus, the system of integral equations (|1.6p can be 
approached by simultaneously minimizing 
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d{yj, Aj x), where d{u,v) is the Kullback-Leibler (KL) 

v{t) 



v{t) In 



u{t) 



v{t) + u{t) 



dt. 



(1.^ 



2 



Throughout this article we will make use of the KL-distance d{v,u) with either u,v £ L^{Q) 
or u, v € L^{^)- 

Remark 1.1 Analog as in (jl.3p . if the assumption A* 1 = 1, j = 0, . . . , N — 1, is satisfied, 
then the first order necessary condition for a minimizer of fj is given by Aj [yj/{Aj x)) = 1, 
and the corresponding multiplicative fixed-point equation reads Pj{x) := x A* {yj/{Aj x)) = x. 

The OS-EM algorithm corresponds to a Kaczmarz type method for solving system (jl.6p 
and can be written in the form 

= P^{xu) = xj ^f'^-^yf^ ds , (1.9) 
Js {AjXk){s) 

where the index < j < N relates to the iteration index k by the formula j = [k] := (k mod 
N). Clearly, the case = 1 corresponds to the standard EM algorithm. 

The cyclic structure of the iteration in (jl.9p is easily recognizable (each cycle consists of 
N steps). Notice that each step within a cycle is an explicit step for solving the fixed point 
equation xA'^)^^ {y[k]/{^[k] ^)) = and can be interpreted as an EM iterative step for solving 
the [fc]-th equation (or block) of system (II. 6p . 

This article is outlined as follows. In Section [2] we formulate a series of assumptions, which 
are necessary for the analytical investigation of the OS-EM method. Moreover, we present 
some basic results concerning the KL-distance. Section [3] contains an analysis of the OS-EM 
iteration (11. 9p . i.e., monotonicity results and consequences concerning the asymptotic behavior 
of the iterations. Section [J] studies the case of noisy data and introduces the loping OS-EM 
method (j4.4p which is a modification of the OS-EM iteration for noisy data. Stability results 
that use discrepancy type principles are stated. In Section [5] we present some numerical 
experiments regarding application of the OS-EM methods to the inversion of the circular 
Radon transform. Section [6] is devoted to final remarks and conclusions. 

2 Assumptions and basic results 

Throughout this article we assume the domains $7 and S in Section [T] to be open bounded 
subsets of R"^, d> 1. The parameter space for investigating system (II. 6p is 

A := {x e L^{n); x>0, x{t)dt = l}, (2.1) 

Jn 

and the starting element xq of iteration (jl.9p is chosen such that xq € A. 

Moreover, we make the following assumptions to the framework introduced in Section [TJ 

(Al) The kernel functions aj : T, x 0, ^ M, j = 0, . . . , N — 1, in (|1.2p satisfy jj, aj{s, t)ds = 1 
for a.e. t € Q; 

(A2) There exist positive constants m and M such that m < aj{s,t) < M a.e. in S x 

(A3) The exact data yj G L^{'^) in (II. 6p satisfy Jj. yj{s) ds = 1; moreover, there exists M' > 
such that yj{s) < M' a.e. in S; 
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(A4) System (jl.7p has a non-negative solution x* G L^(ri), which does not vanish a.e. in Jl; 
moreover, (i(x*,Xo) < oo. 

Assumption (A2) imphes that the operators Aj : L^{Q) — > L^C^) are continuous. More- 
over, any Aj is in L°°(S) and bounded away from zero. This further ensures that ^/ Aj Xk 
has the same properties and then yields that the integrals in ()1.9p are well-defined. 

In the sequel we discuss some basic properties of the KL-distance in (jl.Sp that will be 
needed in the forthcoming sections. This functional plays a key role in the convergence 
analysis of the OS-EM method. For details, we refer the reader to [18^ I17j. 

Lemma 2.1 Let u and v he two functions such that {u,v) is in the domain of the KL- 
distance d{-, •) defined in (jl.Sp . The following assertions hold true: 

(i) d{v, u) > and d{v, u) =0 iff v = u a.e.; 

(ii) \\v - u\\\^ < (§11^11^1 + d{v,u); 
(Hi) The function {v,u) i-^ d{v,u) is convex; 

(iv) Let {vn\ o,nd be given sequences in L^ . If {un} is bounded and lim d{vn,Un) = 0, 

n— >oo 

then lim \\vn — UuWl^ = 0. 

n—foo 

3 The OS-EM method for exact data 

The first result of this section relates to a monotonicity property of the OS-EM iteration. 

Lemma 3.1 Let Assumptions (A1)-(A3) be satisfied, letx G A, and denote Pj{x) = x A*{yj / Aj x). 
Then the following assertions hold true: 

(i) Pj(x) G A and d{Pj{x),x) < f,{x) - fj{Pj{x)), for j = 0, . . . , N - I; 

(ii) If X* A is a minimizer of fj for some < j < N — 1, and d{x*,x) < oo, then 
d{x*,Pj{x)) < oo and fj{x) - fj{x*) < d{x*,x) - d{x*,Pj{x)). 

Proof. Results immediately from |18l Prop. 3.1] applied to the function fj and the corre- 
sponding Pj. □ 

From Lemma [3.1l (i) and Lemma [2.1l (i) we conclude that fj{Pj{x)) < fj{x) and Pj{x) G A. 
Moreover, if x* G A is a solution of (jl.6p with d{x*,x) < oo, then x* minimizes fj, for every 
j = 0, . . . , iV — 1. Lemma |3. II (ii) and the fact that fj{x*) = therefore yield 

fj{x) < d{x*,x)-d{x*,Pj{x)), j = 0,...,iV-l. (3.1) 

In the next lemma we reinterpret the inequalities derived in Lemma 13.11 in terms of the 
OS-EM iteration. 

Lemma 3.2 Let Assumptions (Al)-(A3) he satisfied, and let {xk} be defined hy iteration 
(II. 9|) . Then the following assertions hold true: 

(i) d{xk+i,Xk) < f[k\{xk) - f[k\{xk+i); 
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(ii) If X* € A is a solution of (jl.6p . then f[k]{xk) < d{x*,Xk) — d{x* ,Xk+i)- 

Proof. Results from Lemma l3. II and ()3.ip . □ 

In the next theorem we formulate the main monotonicity results for the OS-EM iteration 
with respect to the KL-distance, as well as convergence results in case the iterations are 
bounded. 

Theorem 3.3 Let Assumptions (Al)-(A3) he satisfied, and the sequence {xk} he defined by 
iteration (11.91). Then we have 



(i) f[k]{xk+i) < f[k]{xk), for every /c G N. 
Moreover, if assumption (A4) is satisfied, then the following assertions hold true: 

(ii) The sequence {d{x* , Xk)} is nonincreasing; 
(Hi) lim fik]{xk) = 0; 

fe— >oo 

(iv) lim d{xk+i,Xk) = 0; 

fc— >oo 

(v) For each < j < N — 1 and p G [1, oo) we have 

lim \\Aj Xj+mN - yj\\LP(T;) = ^ ■ (3-2) 

(vi) If {xf^} is bounded in some U'{Q) space, with p G (l,oo), then it has a subsequence 
which converges weakly in L^{Vt) to a solution of system (|1.6p . 



Proof. Items (i) and (ii) follow from Lemma 13.21 (i) and Lemma |2. II (i). Item (ii) implies 
the existence of > such that lim d{x*,Xk) = /^. Thus, (iii) follows from Lemma 13.21 (ii). 

To prove (iv), notice that (i) and (iii) imply 

lim f[k]{xk+i) = 0. (3.3) 

fc— »oo 

Now, (iv) results from (j3.3p . Item (iii) and Lemma 13.21 (i). 

Next we prove (v). Since fj{x) = d{yj,Aj x), it follows from (iv) that lim d{y[k], Aw-iXk) = 

k — >oo 

0. Consequently, lim d(yj,Aj Xj^mjy) = 0, for every j = 0, . . . , — 1. Now, by applying 

Lemma [2. II (iv) we obtain (j3.2p for p = 1. The case p E (1, oo) follows from [181 Prop. 4.1 and 
Lem. 4.2]. 

The proof of assertion (vi) is divided in several parts: 
(1) Claim: Xk S L°°(J7) for each A: G N. 

Since xq € A by hypothesis, we have m < {AqXq){s) < M a.e. in S by (A2). Consequently, 
it results from (A2) and (A3) that 

^lM|£ifl<^, a.e.inExf^, 

{AoXo){s) m 

and from (jl.9p follows xi E L°°{^}). Part (1) follows by induction if one observes that 
Xk G L~(0) together with (A2) and (A3) imply Xk+i G L°°(0). 
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(2) By hypothesis, the sequence {x^} is bounded in L^'(J7). Therefore, there is a subsequence 
denoted again by {xk}keNj which converges weakly in LP(0) to some z G LP(i}), for some 
pe (l,+oo). 

(3) Conclusion of the proof under a simplifying assumption. 

Let us assume for the moment that, for each fixed < j < — 1, the subsequence {xk}ken 
obtained in Part (2) contains infinitely many indices of the form k = j + mN, m G N. Then, 
for j = 0, we can extract from {x^} a subsequence {x^-} with indices of the form ki = rriiN . 
Obviously, [ki] = for all indices of the subsequence {x^j}, and from (vi) it follows that 
^ki ~^ yo strongly, and thus weakly in LP(5]). Since is continuous from LP{Vl) to L'P(S) 
due to (A2), it is weakly continuous. Part (2) implies that A^Xk^ Aq z weakly in LP(S). 
From the uniqueness of weak limits it follows that •2 = 2/0- 

By repeating the argumentation for j = 1, . . . , N — 1, we conclude that Aj z = yj, for j = 
1, . . . , N — 1, thus proving that z is a solution of system ()1.6p . 

(4) Conclusion of the proof in the general case. 

If the assumption in part (3) does not hold, then there must be at least one < jo < N — 1 
such that the subsequence {xk}keN obtained in part (3) contains infinitely many indices of 
the form k = Jq + mN, m G N. Arguing as in part (4) we obtain a subsequence {xk^}, with 
indices of the form ki = jo + rriiN, i € N, such that Aj^^Xk^ yjo weakly in L^(S), and 
conclude that Aj^ z = yj^ . 

Now, let us consider the subsequence {xk with indices of the form ki + 1 = (jo + 1) + mjA^, 
Item (iv) implies that d{xki-\-i, Xk^) ^ 0, as i — > oo. Since both subsequences 
{xfc._|_i}j, {xki}i are in L'p{Q) and bounded (part (1) above), it follows from Lemma l2.ll (iv) 
that — XkiWhP — 5- as fcj ^ oo. Therefore, x^.+i — > z weakly in and from the 

continuity of follows that .4yg_|_i] x^.+i — > .4[jg_|_i] z weakly in L^'(S). Moreover, (v) 

yields Xk^+i — > Vjo+i weakly in WiYi). Then we conclude that ^[jg+i] z = yy^+i]. 

Repeating the argumentation for the subsequences {xfc-_|_2}, • • . , {xk^+N-i}-, we conclude that 
Aj z = yj, for every j, proving that z is a solution of system (11. 6|) . □ 

Remark 3.4 We can interpret Theorem \3.3\ (v) as follows: If we consider the subsequence 
{xj+mAf}meN formed by the j-th component of each cycle of the OS-EM iteration (where 
< J ^ N — 1), then the L^-norm of the residual corresponding to this subsequence converges 
to zero. 

Moreover, Theorem \ 3.^ (iv) guarantees that, given any two "consecutive" subsequences 
{xj+mN}m(in and {x(j+i)+mAf }meN, we have 

lim d(x(j+i)+mAfi a;j+mAr) = 0, 

m— >oo \j / • 

for each j = {),... ,N — 2. 

4 The loping OS-EM method for noisy data 

Our next goal is to modify the OS-EM iteration by introducing a relaxation parameter, and to 
investigate monotonicity and stability results for this modified iteration (the so called loping 

'^Notice that {ajfc^+i} may contain elements which do not belong to the convergent subsequence {a;*;} obtained 
in part (3), while {xk^} is a subsequence extracted from {xfc}. 
^Notice that Xk- z weakly in 
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OS-EM method) in the case of noisy data. As remarked in [8], "With noisy data though, 
inconsistent appHcations (of discrete OS-EM - authors' note) result." 

We aim at characterizing the loping OS-EM method as an iterative regularization method 
in the sense of [3]. 

For the rest of this section we assume that the right hand side of (11. 6p is not exactly 
known. Instead, we have only approximate measured data 7/| E satisfying 



\\yj-y%^ < S,, j = o,...,N-i. (4.1) 

We denote 5 := {6o, ...,6]\f-i)- 

In this noisy data case we are interested in finding an approximate solution for the system 

Ajx = yj, j = 0,...,N-l. (4.2) 

The following assumptions are required for the analysis: 

(A5) The noisy data G L^{'^) satisfies Vji^) '^^ = 1- 

(A6) There exist Mi, mi > such that Mi > y| > mi a.e. in S. 

Also necessary for the analysis are the following functions associated to the equations of 
system (|42]) 

//(x) := ljy'^{s)\nj^^^-y'^is) + {Ajx){s)]ds. (4.3) 

Notice that f^{x) = d{y^j,Aj x). 

The loping OS-EM iteration for the inverse problem ()4.2p with noisy data is defined by 

xt+i = xiuJk (4.4a) 

where 



f^'-T^"' 4l(-^)>-^%l. (4.4b) 



1 , else 
The constants r and 7 in (|4.4bp are chosen such that 



r > 1, 7 =max<^ In— , In \\, (4.5) 

where m, M, mi. Mi are the positive constants defined in (A2) and (A6). 

Remark 4.1 It is worth noticing that, for noisy data, the iteration in (|4.4p is much different 
from the iteration in (jl.Op ." The relaxation parameter loi^ effects that the iterates defined in 
()4.4ap become stationary if all components of the residual vector d{y^j^Y A]j.^x^^ fall below a 
pre-specified threshold. 

Another consequence of using these relaxation parameters is the fact that, after a large 
number of iterations, Wfc = 1 for some k within each iteration cycle. Therefore, the computa- 
tional evaluation of the adjoint operator 
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might be loped, making the loping OS-EM iteration in (j4.4p a fast alternative to the OS-EM 
method. 

In the case of noise free data, i.e. 5j = in (14. ip . we choose uj^ = P^^{x^j^ = P]j.-^{xk) and 
the loping OS-EM iteration (Oj) reduces to the OS-EM method (fL9]) . 



In the sequel we prove a monotonicity result for the loping OS-EM iteration in the case 
of noisy data. First however, we derive an auxiliary estimate. 

Lemma 4.2 Let assumptions (A1)-(A5) hold true. Moreover, let i/j, Sj be given as in (14. ip . 
with > for some < jo < — 1. Then we have 

/[l](4)-%[fc],4]) ^ d{x\xi)-d{x\xi^,). (4.6) 
for a// A; € N with [k] = jo . 

Proof. Since (Al)-(A3) are satisfied, we argue as in the proof of [181 Prop. 5.2] to conclude 
that for every v, w (z A, and < j < — 1 the inequality 

d{Pj{w),Pl{v)) < d{yj,y'j)+d{Pj{w),v)-d{Pj{w),w)+fj{w)-fj{v), 

holds true. Therefore, given A; S N with [k] = jo, ()4.6p follows by taking j = [k], w = x*, 
V = x^f,, and by observing that {x* ) = x* . □ 

Proposition 4.3 Let assumptions (Al)-(A6) hold true andr, 7 be defined as in (j4.5p . More- 
over, let yj, 6j be given as in (14. ip with 6j > for j = 0, . . . , N — 1. Then the sequence {x|} 
defined by iteration ()4.4p satisfies 

d{x\xi^^) < d{x*,xi), ken. (4.7) 

Proof. If f^f^-^ixf.) < TjS^k], then u;^ = 1 by ()4.4bp . Therefore, xf._^-j^ = xf. and (|4.7p follows 
with equality. If f^j^jix^) > T^S^k]i notice that a simple calculation yields 



d{x*,xi)-d{x*,xi+,) > /[i](4) + JJy[k]{s)-yf,]{s)]ln(^ 
from (j4.6p . Therefore, ()4.7p follows from 

4] (4)+ / [y[k]-yfk]]^n{^)ds > 



ds 



> 



> f[k]i4) - ^k] max {1 In ^1 , I In 

>/[l](4)-7'^[fc] (4.8) 

> (r - l)7 5[fc] 
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together with (|4.5p . To obtain the inequahties above we used (|4.ip . (j4.5p . (A2) and (A6). □ 

Proposition 14.31 gives us a hint on how to choose the stopping rule for the loping OS-EM 
iteration. That is, we stop the iteration at 

kl := mm{mN G N ; x^^ = x^^+i = • • • = x^^+^.J . (4.9) 

In other words, fcf is the smallest integer multiple of N such that 

Xks = Xks+i = ■■■ = XkS+N-i ■ (4-10) 

In the sequel, we prove that the stopping index in (14. 9p is well defined and that the 
corresponding iterations stably converge to a solution of the system, if they are bounded in 
some U' space with p € (1, +oo). 

Theorem 4.4 Let assumptions (A1)-(A6) he satisfied, and G N 6e chosen according to 
(j4.9p . Then the following assertions hold true: 

(i) The stopping index k^ defined in (j4.9p is finite; 
(a) More precisely, kl = 0{5^^^), where Jmin := min{(5o, . . . , ^at-i}; 
(Hi) d{y^j,Aj xls) < T-f6j, for every j = 0, . . . , iV - 1; 
(iv) For every p € [1, +oo) and every j = 0, . . . , N — 1 we have 

limJAjxls -yillLP{E) =0. 

(v) Let {d'- := ((5q, . . . , (5^_j^)}/gN be a sequence in (0, oo)^ with lim M = 0, for each < 
j < N — 1. Moreover, let {y' := (i/q, . . . , y5v-i)}'eN be a sequence of noisy data satisfying 

llyj-yjilLi < 4' i = o,...,A^-i, /GN, (4.11) 

and kl := kf = ktf{6\y'') denote the corresponding stopping index defined in ()4.9p . If 
the sequence {x^y is hounded in some LP{Q) space, with p E (1, +oo), then it has a 
suhsequence which converges weakly in LP{0,) to a solution of system (|1.6p . 

Proof, (i) Assume by contradiction that k^ is not finite. Then it results from (j4.9p that 
^k+i 7^ -^fe least once in each cycle of iteration (j4.4p . Hence for every m S N there exits 
jm e {0, . . . , iV - 1} such that 

fjMm+mN) > rjS,^ ■ (4.12) 
From (j4.8p in the proof of Proposition 14.31 it follows 

dix*,xi)-d{x*,xi+,) > max{4](xf)-7<^[fc],0}, kGN. 
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Summing up this inequality for k = 0, . . . , IN — 1 implies H 

IN-l 

d{x*,xo) - d{x*,xf^) > ^ max{/j'^](xfc) - 75[fc], 0} 

fe=o 

I N-l 
m=0 j=0 

Then, it follows from (|in^ 

z z 

m=0 m=0 

(4.13) 

However, due to (|4.5p . the right hand side of (j4.13p becomes unbounded as Z — > oo, contra- 
dicting (A4). Therefore, kl must be finite. To prove (ii), it is enough to take / = ki/N e N 
in (|4l3]) and obtain k^ < N d{x* , xq) / {{t - l)7(5min)- 
To prove (iii), we assume by contradiction that 

for some < jo < iV - 1. Thus, it results from (liTOl) that //^(x^^^^J > T^iSj^. Therefore, it 
follows from ()4.8p in the proof of Proposition 14.31 that 

= d{x\xls^^^)-d{x\xis^^^^,) > f^,{xis^^J-j6,, > {T-lh6,,, 

which contradicts (j4.5p . 

(iv) and (v) The proofs follow the lines of the proof of Theorem 13.31 (v), (vi). □ 

Remark 4.5 (Stability for noisy data in L^(S)) When dealing with inverse problems, bounds 
for the noisy data are most commonly given in the L^-norm, i.e. the approximate measured 
data Hj S L'^i'^) is assumed to satisfy 

\\yj-y%^ < Sj, j = o,...,N-i, (4.14) 

instead of (14. ip . In this case, the loping OS-EM iteration is defined by ()4.4p . where the "loping 
condition" f^i,j{xf.) > Tjd^i^j in ()4.4bp is substituted by 

f[k](.4) > r5[k]\\Hytk]/iAk]4mL^- (4.15) 

Under this assumptions it is possible to state a stability result, similar to the one in Theorem 
\4.4\ (iv). One argues as follows: 

• First of all, notice that monotonicity of the error with respect to the KL-distance (as in 
()4.7p ) follows when using the Cauchy-Schwarz inequality in i^^(S) to derive the estimate 
( compare with (j4.8p ) 

f[k]i4) + j hk] - y[k]\ HyU/iAk] 4)) ds>{T- 1) || Hyl^/iA^u] 4))\\l^ ■ 

^Notice that Xq = xo- 
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• By defining the stopping index as in ()4.9p , its finiteness can be proven analogously as 
in Theorem \4.4\ (i)- Moreover, the following estimate holds true (compare with Item (Hi) 
of Theorem \4-4^ 

d{y^j,Aj xls) < t6j\\ ln{yy{A[k] j = 0, . . . , - 1 . (4.16) 

• In Theorem \4-4\ ('^'V), if one substitutes the assumption (j4.1ip by \\yj — yjHi^ < ^j, 
J = 0, . . . , — 1, / G N, then the proof of the stability result carries on with analogous 
argumentation. 

Notice that the estimate in (j4.16p allows for the following interpretation: The loping OS- 
EM iteration should be stopped at the index ( an integer multiple of N) when for the first 
time (|4.16p is satisfied within a whole cycle. 

The advantage of using this stopping rule resides on the fact that no quantitative in- 
formation on the constants m, M, mi, Mi is required to compute the iteration. In other 
words, the constant 7 is not required neither to test the "loping condition" (|4.15p nor to verify 
the stopping rule based on (j4.16p . This is obviously not the case if the "loping condition" 
/[fcl(^fc) ^ ''"7^[fc] (!4.4bp is to be implemented. 



5 Numerical example 

In this section we compare the numerical performance of our loping OS-EM method with the 
OS-EM and EM methods. As benchmark problem we use a system of linear equations for the 
circular Radon transform. The inversion of the circular Radon is relevant for the emerging 
photoacoustic computed tomography [12l [151 EQl [22] . 

Let e < 1 be some small positive number, let i7 := i?i_e(0) C denote the disc with 
radius 1 — e centered at the origin, set 

and let $ : M ^ M be a continuous nonnegative function with supp(<I>) = [— e, e] and jjg <I> = 1. 
Our aim is the stable solution of (jl.6p . with Aj x := ^ *r {Mj x), where 

rN f 

{M-j x){ip,r) := / x{{cos ip, sin ip) -\- ru>) du , (99, r) € S,- , (5.1) 

is the circular Radon transform restricted to , and ^*ry = T^ y denotes the convolution of 
and y . In (j5.ip . x is considered as an element in L^(E?) by extending it with zero outside 
of 17. 

One verifies that the operators Aj can be written in the form (jl.7p . with s = (ip, r) and 
aj{t, r) = $ (|(cos ip, sin 99) — t| — r) , j = 0, . . . , — 1 . 
Moreover, the adjoint of Aj is given by ^* y = Bj{(^ *r y), where 

{Bjy){t) = —j y{\t - {cos p,sm.ip)\) dp , (5.2) 

J2jn/N 

is the circular backprojection. Hence ^* 1 = 1 and the operators Aj satisfy assumption (Al). 
However, since aj are not bounded from below, Aj do not satisfy (A2). 
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Remark 5.1 For any positive A, the operators 



Af^ X := \— {AiX + \ I x] , j = 0,...,N -1, 



l + A|Sj| V Jn 

clearly satisfy (A2). Since {Af^fy = {A*y + \j^,y)/{l + X\T.j\) we have {Af'')*l = I, 
proving that (Al) is also satisfied. 

Therefore, we shall consider for the rest of this section the system of equations 

Afx = yf:=-^^^^^(^, + \j^ y,^ , j = 0, . . . , iV - 1 . (5.3) 

The identity JqX = J^x A* 1 = jj, AjX implies that x is a solution of ()5.3p if and only if x 
satisfies Aj x = yj. 



If noisy data yj with \\yj — VjWl'^ < Sj are available, then 



where yj^^'^ is defined in the same way as yj^\ with yj replaced by Therefore, the loping 
OS-EM iteration with noisy data applied to system ()5.3p reads as 



1+A|S,| U*A^w4.+aJ ' 1+A >^lO[k], V ) 

1 , else . 



Here we made use of the fact that the initial guess satisfies = 1, which implies 

ku, Ak] = 1 for every k. 



Remark 5.2 Iteration \5.4^ assumes continuous data yj G L^{Tij), whereas in practical ap- 
plications only discrete data are available. In the following we assume that data 

y%,ir] :=yj(<p[v],r[v]), G . . . , {j + 1)N^ - 1} x {0,...,iV,}, 

are given, with <p[iip] '■= 2i^TT/N^, r[ir] := 2ir/Nr, and N^p, Nr + 1 denoting the number of 
samples of y^ in the angular and radial variable, respectively. 

In the numerical implementation Mj, Bj, and d are replaced (as described below) with 
finite dimensional approximations Mj, Bj, 1$, d, and (j5.4p is approximated by 

xi+i(t[i]) ~ H := xi[i\ui,[i\ , i G {0, . . . ,iVt}2 

r Bwl,.+A / <;+A X d(4j+A,I,.Mff4+A) (^) 
Uk := < l+47rA/iV l,I<i,M[fc] xi+\J ' 1+A ^ ^'"[k] ' 

[ 1 , else . 

Here = (yj[v, v])v,i,, = (a;|[i])i, ^fc = (i^fcH)i, and = -(1,1) + 2i/Nt with 
(A'f + 1)^ denoting the number of samples in the variable t. 
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1. The discretized circular Radon transform 

is obtained by replacing x in (j5.ip with the bilinear spline T{x) satisfying T{x){t[i]) = 
x[i], and approximating the resulting integrals over the with the trapezoidal rule. 
This leads to 

2 3r[ir]Nt 

iMjx)[i^,ir] = — T{x)(^a[i^]+r[ir]u[i^]), (5.6) 

where <T[i<^] := {cos (p[i^], sin Lp[i^]), u[iu)] := {cos{2m^ / Nt), sm{2m^ /Nt)), and 3r[ir]Nt 
is the number of supporting points when applying the trapezoidal rule. 

2. Assuming that e = 2K/Nr for some iiT G N, the convolution Z$ y = <I> *f. y is approxi- 
mated by 

(I$2/)[V,V] = — ^{'^ii-ir)/Nr)y[ip,i'r], 

where y[iy,,ir] '■= for ir outside {0, . . . , Nr}. 

3. The discretized back-projection Bj : R^^xi^r+i) ^iNt+i)x{Nt+i) defined by 

^ (i+i)A^^-i 
(Bjy)[i]:=— Yl Tr{y){i^,m-cT[i^]\), 

if t[i] G 0, and setting (Bjy)[i] := for t[i] Q. Here Tr{y) denotes the piecewise 
linear spline in the second variable satisfying Tr{y){i^,r[ir]) = y[iip,ir]- 

4. Finally, the discrete KL-distance is defined by 

4^ (j+i)Af^-i Nr v[i i] 

for -j;,it G RA^v=x(A^r+i)_ 

Remark 5.3 (Numerical Complexity) Assuming Nt = Nr, the numerical complexity for 
performing one iteration cycle (which consists of N subsequent steps in i5. 5\) ) is ©(A'^angic-^*^)- 
Here A'^angic = N^N corresponds to the overall angular data samples, which is independent 
of N in practice. Therefore in the following we always compare the reconstruction error in 
dependence of the number of iteration cycles. 

In the following numerical examples we apply the (loping) OS-EM iteration with = 1 
(corresponding to the EM algorithm), A = 5, = 10, and A = 20 subsets. The original 
phantom x* (the exact nonnegative solution) is shown in the left picture in Figure [1] and 
consists of a superposition of characteristic functions. Note that a similar phantom was 
reconstructed in [8] where the OS-EM technique was introduced. The data y,-, shown in the 
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Figure 1: Original phantom x* (left) and simulated data {AjX*)j. 



right picture in Figure[Tl were calculated numerically by (j5.6p for A'^angic = N^pN = 100 angular 
samples. In order to avoid inverse crimes, much larger Nt is used for the data simulation as 
for the application of the loping OS-EM iteration. In all examples a;o = a;Q = l/((l — e)^7r) is 
used as initial guess and the parameters e and A are chosen to be 0.02 and 0.01, respectively. 

The iterations of the OS-EM method applied to exact data with Nt = = 100 and 
different values of N are depicted in Figure [2j It can be seen that the 5-th iterate with EM 
has similar quality as the 1-th iterate with OS-EM for N = 5. As a rough rule one can say 
that making cycles with the EM algorithm leads to an improvement similar to 1 cycle with 
the OS-EM algorithm. This can also be recognized in the left image in Figure [3l where the 
evolution of the error is depicted with respect to the KL-distance. 

In order to investigate the dependence of the OS-EM iteration on the discretization level, 
we repeated the experiment with Nt = N,. = 200. The right image in Figure [3] shows the 
corresponding logarithmic error. As expected, the error is relatively independent on the 
discretization. 

In the case of noisy data we apply the loping OS-EM iteration ()5.5p . The noisy data 
is created by adding 5% Poisson distributed noise to the simulated data yj, such that 
47T/{NrN^)Y:\yjK,ir] " yj[v, v]| = 0.05. 

Remark 5.4 Our numerical experiments show that, for large 6 and r = 1, far too many 
iterations are loped. A significant improvement can be obtained if t = t{5) is chosen in 
dependence of the noise level, with t{5) < 1 for large 5 and t{5) converging to some Too > 1 
for (5 — > 0. It is clear that the asymptotic convergence analysis (for 6^0) remains valid in 
such a situation. 

The reconstruction for noisy data with Nt = Nr = 100 are depicted in Figure HI For com- 
parison purposes, results of the OS-EM iteration (without loping strategy) are also included. 
The loping OS-EM is automatically stopped according to (14. 9p whereas their non-loping coun- 
terparts are stopped after the cycle with minimal error d{x*,xf.), which is not available in 



14 




Figure 2: Exact data experiment: Iterates for Nt = = 100 with = 1 (left), = 5 
(middle) and = 10 (right) after 1, 5 and 25 cycles. 

practice. All reconstructions are quite comparable. Figure [5] shows the evolution of the er- 
ror with respect to the KL-distance. In this figure one also notices the semi-convergence of 
the non-loping iterations, which happens typically when applying non-regularized iterative 
schemes to ill-posed problems [2| [71 [TO]. 

An inspection of Figure [5] shows that the regularized solution of the loping OS-EM methods 
(automatically stopped) have errors comparable to the optimal solution of their non-loping 
counterparts when stopped after the cycle with minimal error (which is not available in the 
practice). Figure [6] shows the number of actually performed iterations. Table [T] summa- 
rizes run times and errors with Nt = = 100, A^angle = 100 (with non-optimized Matlab 
implementation on HP Notebook with 2 GHz Intel Core Duo processor). 
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log d{x*,Xk) for Nt = 100 log d{x*,Xk) for Nt = 200 




5 10 15 20 25 '^'^0 5 10 15 20 25 



Figure 3: Exact data experiment: Logarithmic plots of iteration error with respect to the 
Kullback-Leibler distance for Nt = Nr = 100 (left) and Nt = Nr = 200 (right). 





N 


^cycl 


time (sec) 


d{x*,xi) 


loping OS-EM 


10 


4 


13.4 


0.022 


OS-EM 


10 


3 


9.2 


0.022 


loping OS-EM 


20 


4 


13.4 


0.024 


OS-EM 


20 


2 


6.3 


0.024 



Table 1: Comparison of the performance of different iterative methods. The non- loping 
iterations are stopped after the cycle with minimal error, whereas the loping OS-EM are 
automatically stopped according to (j4.9p . 

6 Conclusions 

This article is devoted to the investigation of OS-EM type algorithms for solving systems 
of linear ill-posed equations. We focus on showing regularization properties of the proposed 
methods. 

In the case of exact data, our approach originates an algorithm analog to the OS-EM 
iteration. We are able to prove monotonicity results with respect to the Kullback-Leibler 
distance as well as weak convergence in case of boundedness of the iterations. In the noisy 
data case, we propose a loping OS-EM iteration which differs from the OS-EM method due 
to the introduction of a loping strategy. This loping strategy renders the proposed iteration a 
regularization method. We prove monotonicity of the iterates and study stability properties 
of our method. 

What concerns numerical effort, we conjecture that the loping OS-EM algorithm is at 
least as efficient as the well established OS-EM method. The numerical experiments with 
(15. 5|) for inverting the circular Radon transform support this conjecture. In the case of exact 
data, (j5.5p reduces to a discretized version of the continuous OS-EM iteration applied to the 
system (|5.3p . However it is slightly different to the discrete OS-EM iteration of [8] since Bj is 
not the exact transpose of Mj. Moreover, opposed to [8] our continuous convergence analysis 
applies independent on the discretization level. 
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Figure 4: Noisy data experiment: Iterates without loping (top line) and with loping (bottom 
line). The loping iterations are stopped automatically whereas their non-loping counterparts 
are stopped at the iteration cycle where d{x*,xf.) is minimal (which is not available in prac- 
tice) . 
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Figure 5: Noisy data experiment: Evolution of the relative error log d(a;*, a;|) for loping and 
non-loping OS-EM iterations. 




References 

[1] A. De Cesaro, M. Haltmeier, A. Leitao, and O. Scherzer. On Steepest-Descent-Kaczmarz 
methods for regularizing systems of nonlinear ill-posed equations. Appl. Math. Comput., 
202:596-607, 2008. 

[2] P. Deuflhard, H. W. Engl, and O. Scherzer. A convergence analysis of iterative meth- 
ods for the solution of nonlinear ill-posed problems under affinely invariant conditions. 
Inverse ProbL, 14:1081-1106, 1998. 

[3] P. P. B. Eggermont and V. N. LaRiccia. Maximum penalized likelihood estimation and 
smoothed EM algorithms for positive integral equations of the first kind. Numer. Fund. 
Anal. Optim., 17(7-8) :737-754, 1996. 

[3] H. W. Engl, M. Hanke, and A. Neubauer. Regularization of inverse problems, volume 
375 of Mathematics and its Applications. Kluwer Academic Publishers Group, Dordrecht, 
1996. 



18 



[5] M. Haltmeier, R. Kowar, A. Leitao, and O. Scherzer. Kaczmarz methods for regularizing 
nonlinear ill-posed equations II: Applications. Inverse Prohl. Imaging, 1:507-523, 2007. 

[6] M. Haltmeier, A. Leitao, and O. Scherzer. Kaczmarz methods for regularizing nonlinear 
ill-posed equations I: convergence analysis. Inverse Prohl. Imaging, 1:289-298, 2007. 

[7] M. Hanke, A. Neubauer, and O. Scherzer. A convergence analysis of Landweber iteration 
for nonlinear ill-posed problems. Numer. Math., 72:21-37, 1995. 

[8] H. M. Hudson and R. S. Larkin. Accelerated image reconstruction using ordered subsets 
projection data. IEEE Trans. Med. Imag., 13:601-609, 1994. 

[9] Alfredo N. lusem. A short convergence proof of the EM algorithm for a specific Poisson 
model. Rebrape, 6(l):57-67, 1992. 

[10] B. Kaltenbacher, A. Neubauer, and O. Scherzer. Iterative Regularization Methods for 
Nonlinear Ill-Posed Problems, volume 6 of Radon Series on Computational and Applied 
Mathematics, de Gruyter, Berlin, 2008. 

[11] R. Kowar and O. Scherzer. Convergence analysis of a Landweber-Kaczmarz method 
for solving nonlinear ill-posed problems. /// posed and inverse problems (book series), 
23:69-90, 2002. 

[12] P. Kuchment and L. A. Kunyansky. Mathematics of thermoacoustic and photoacoustic 
tomography. European J. Appl. Math., 19:191-224, 2008. 

[13] H. N. Miilthei and B. Schorr. On an iterative method for a class of integral equations of 
the first kind. Math. Methods Appl. Sci., 9(2): 137-168, 1987. 

[14] H. N. Miilthei and B. Schorr. On properties of the iterative maximum likelihood recon- 
struction method. Math. Methods Appl. Sci., ll(3):331-342, 1989. 

[15] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. Experimental evaluation of 
reconstruction algorithms for limited view photoacoustic tomography with line detectors. 
Inverse Probl, 23(6):81-94, 2007. 

[16] R. A. Redner and H. F. Walker. Mixture densities, maximum likelihood and the EM 
algorithm. SI AM Rev., 26(2):195-239, 1984. 

[17] E. Resmerita and R. S. Anderssen. Joint additive Kullback-Leibler residual minimization 
and regularization for linear inverse problems. Math. Methods Appl. Sci., 30(13):1527- 
1544, 2007. 

[18] E. Resmerita, H. W. Engl, and A. N. lusem. The expectation-maximization algorithm for 
ill-posed integral equations: a convergence analysis. Inverse Problems, 23(6):2575-2588, 
2007. 

[19] W. H. Richardson. Bayesian-based iterative method of image restoration. J. Opt. Soc. 
Am., 62:55-59, 1972. 



19 



[20] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen. Variational 
Methods in Imaging, volume 167 of Applied Mathematical Sciences. Springer, New York, 
2008. 

[21] Y. Vardi, L. A. Shepp, and L. Kaufman. A statistical model for positron emission 
tomography. J. Amer. Statist. Assoc., 80(389) :8-37, 1985. With discussion. 

[22] M. Xu and L. V. Wang. Photoacoustic imaging in biomedicine. Rev. Sci. Instruments, 
77(4):041101, 2006. 



20 



